%% 0.1 一维原子链模拟程序OneDimenChain
%
% 处理模型：弹簧-质量小球模型~o~~o~~o~~o~
% 采用新单元模型：~o~
%
% 积分格式：【速度Verlet方法】
% 作者：hellobei
% 日期：2020-03-23
%
% 矩阵全部采用稀疏矩阵存储
%
% 用法：
%   1. 首先：运行ODCpre.m设置参数。
%   2. 动画显示结果：运行ODCpost.m
%   2. 初始条件在4中设置
%   3. 边界条件在5中设置。

%% 1. 准备工作
tic     % 开始计时

fprintf('[ODCmain] Begin to calculate.\n');
fprintf(fid,'\n[ODCmain] Begin to calculate.\n\n');

%% 2. 生成质量阵M, 刚度阵K, 阻尼阵C
[M,K,C,m,k,nBall,DOF]=func_model(L,ctn,cm,ck);

fprintf(fid,'[ODCmain] DOF = %d\n',DOF);
fprintf(fid,'[ODCmain] m = %0.3f\tk = %0.3f\n',m,k);

%% 3. 分析步step
iframe=1; %frame编号（从1开始）
nframe=Tt/dt+1; %总的frame帧数-1

fprintf(fid,'[ODCmain] Frame number nt is %d\n',nframe);

phyTime=0; % phyTime=(iframe-1)*dt;
iout=1; %常变量存储矩阵A的列编号（从1开始）

%% 4. 初始条件（手动设置：边界条件若存在，则在此修改第一部分 *）
[Q,a,da]=func_iniCondition(DOF); %设置初始条件

a(1)=func_load(iframe,w,dt); %边界条件（手动设置*）

dda=M\(Q-K*a-C*da); %加速度自动确定

outA(:,iout)=a; %为画图方便，保存位移变化的时间历程，共计n+1列

%% 5. 能量计算
[se,ke]=func_energy(m,k,nBall,a,da);% 计算此时刻的能量

[totalE,doneE,doingE,doneP,doingP]=func_recordEnergy_ini(Tt); %记录能量

%% 6. Velocity-Verlet格式求解（手动设置：边界条件若存在，则在此修改第二部分 *）
%% 输出信息
fprintf(fid,'[ODCmain] Integrator Option is %d\n',option);

[c0,c1,c2,c3,c4,c5,c6,c7,interMatrix] = func_Inte_ini(option,dt,M,K,C);

% 辅助参数，为了输出方便
count=0;
fprintf('[ODCmain] Now CPU time is ');

while iframe<nframe
    %% 进入下一步
    iframe=iframe+1;
    
    %% 积分求解器 核心
    [a,da,dda] = func_Inte(option,c0,c1,c2,c3,c4,c5,c6,c7,dt,a,da,dda,Q,M,K,C,interMatrix);
    
    %% 边界条件（手动设置*）
    if((iframe-1)*dt<=t_cut)
        a(1)=func_load(iframe,w,dt);
    end
    
    %% 场变量输出
    if (rem(iframe-1, dframe)==0)
        iout=iout+1;
        outA(:,iout)=a; % 位移场
    end
    
    %% 能量
    if (rem((iframe-1)*dt, 1)==0)
        [se,ke]=func_energy(m,k,nBall,a,da);
        
        phyTime=phyTime+1; % ClassicLength=WaveSpeed*phyTime
        
        [totalE,doneE,doingE,doneP,doingP]=func_recordEnergy(se,ke,phyTime,t_cut,ctn,totalE,doneE,doingE,doneP,doingP);
    end
    
    % 输出计算时间
    costtime=toc;
    fprintf(repmat('\b',1,count));
    count=fprintf('%0.2fs, [%6.2f%%], [%6d/%6d]',costtime,(iframe/nframe)*100, iframe,nframe);
end

% 保存工作区
save(workname_path)
fprintf(fid,'\n[ODCmain] Save to %s.mat\n',workname_path);

% 结束
fprintf(fid,'\n[ODCmain] Calculate over. Cost %0.2fs\n',costtime);
fprintf(fid,'\n=================================================\n');
fprintf('\n[ODCmain] Calculate over. Please run ODCpost.m\n');
